function [FEDec_Y,FEDec_X,Std_Y,Std_X]=ForecastErrorDecomposition(Mat_A,Mat_B,Mat_C,MaxT)
% X_{t+1}=A*X_{t}+B*Shock_{t+1}, Y_{t}=C*X_{t}
Num_X       =   size(Mat_A,1);
Num_S       =   size(Mat_B,2);
Num_Y       =   size(Mat_C,1);

FEDec_X     =   zeros(Num_X,Num_S,MaxT);
FEDec_Y     =   zeros(Num_Y,Num_S,MaxT);
DecVar_X    =   zeros(Num_X,Num_S,MaxT);
DecVar_Y    =   zeros(Num_Y,Num_S,MaxT);
Std_X       =   zeros(Num_X,MaxT);
Std_Y       =   zeros(Num_Y,MaxT);

% For Dt=1, or the impact of eps_{t+1} on e_hat(X_{t+1}|X_{t})
TempCoefMat_X           =   Mat_B;
TempVar_X               =   TempCoefMat_X.^2;
TempCoefMat_Y           =   Mat_C*TempCoefMat_X;
TempVar_Y               =   TempCoefMat_Y.^2;
DecVar_X(:,:,1)         =   TempVar_X;
DecVar_Y(:,:,1)         =   TempVar_Y;
Std_X(:,1)              =   sum(DecVar_X(:,:,1),2);
FEDec_X(:,:,1)          =   bsxfun(@rdivide,DecVar_X(:,:,1),Std_X(:,1));
Std_Y(:,1)              =   sum(DecVar_Y(:,:,1),2);
FEDec_Y(:,:,1)          =   bsxfun(@rdivide,DecVar_Y(:,:,1),Std_Y(:,1));
% For Dt>1
if MaxT>=2
    for tt=1:MaxT-1
        TempCoefMat_X           =   Mat_A*TempCoefMat_X;
        TempVar_X               =   TempCoefMat_X.^2;
        TempCoefMat_Y           =   Mat_C*TempCoefMat_X;
        TempVar_Y               =   TempCoefMat_Y.^2;
        DecVar_X(:,:,tt+1)      =   TempVar_X+DecVar_X(:,:,tt);
        DecVar_Y(:,:,tt+1)      =   TempVar_Y+DecVar_Y(:,:,tt);
        Std_X(:,tt+1)           =   sum(DecVar_X(:,:,tt+1),2);
        FEDec_X(:,:,tt+1)       =   bsxfun(@rdivide,DecVar_X(:,:,tt+1),Std_X(:,tt+1));
        Std_Y(:,tt+1)           =   sum(DecVar_Y(:,:,tt+1),2);
        FEDec_Y(:,:,tt+1)       =   bsxfun(@rdivide,DecVar_Y(:,:,tt+1),Std_Y(:,tt+1));
    end
end
Std_X       =   sqrt(Std_X);
Std_Y       =   sqrt(Std_Y);